function D=osl_forward_model(S)

% D=osl_forward_model(S)
% 
% Required inputs:
% S.D: SPM MEG object filename
% S.sMRI: structural MRI nii file name (set S.sMRI=[] or '' to use template
% structural)
% S.useheadshape: set to 0 or 1 to indicated if the headshape points should
% be used in the registration
%
% Check the output using:
% spm_eeg_inv_checkmeshes(D);
% spm_eeg_inv_checkdatareg(D);
% spm_eeg_inv_checkforward(D, val);
%
% MW

global OSLDIR;

try, useheadshape=S.useheadshape; catch, error('S.useheadshape not specified'); end;
try, spmfilename=S.D; catch, error('S.D not specified'); end;
try, sMRI=S.sMRI; catch, error('S.sMRI not specified'); end;
try, do_neuromag_grad_baseline_correction=S.do_neuromag_grad_baseline_correction; catch do_neuromag_grad_baseline_correction=1; end;
try, correct_grads=S.correct_grads; catch, correct_grads=1; end;
    
val=1;
do_eeg=0;

if(strcmp(sMRI,''))
    sMRI=[OSLDIR '/std_masks/MNI152_T1_1mm.nii'];
    
    warning(['Using ' sMRI ' as structural']);
end;

D = spm_eeg_load(spmfilename);

if sum(ismember(unique(D.chantype),'MEGGRAD')) && do_neuromag_grad_baseline_correction % HL Mod 1.3 - Prevents correction of CTF data.
    warning('Axial Gradiometers Detected! Are you using CTF data, If so you should not be correction gradiometers. Setting do_neuromag_grad_baseline_correction to 0.')
    do_neuromag_grad_baseline_correction=0;
end

[ok, D] = check(D, 'sensfid');

[D,val] = spm_eeg_inv_check(D, val);
val=1;
if ~isfield(D, 'inv') || ~isfield(D.inv{val}, 'comment')
    D.inv = {struct('mesh', [])};
    D.inv{val}.date    = strvcat(date,datestr(now,15));
    D.inv{val}.comment = {''};
else
    tinv = struct('mesh', []);
    tinv.comment = D.inv{val}.comment;
    tinv.date    = D.inv{val}.date;
    D.inv{val} = tinv;
end

D.inv{val}.mesh = spm_eeg_inv_mesh(sMRI, 2);


%-Check meshes and display
%spm_eeg_inv_checkmeshes(D);

%D.save;

%%%%%%%%%%% coregister/datareg

%D = spm_eeg_load(spmfilename);

meegfid = D.fiducials;
mrifid = D.inv{val}.mesh.fid;

meeglbl = meegfid.fid.label;
mrilbl = mrifid.fid.label;

newmrifid = mrifid;
newmrifid.fid.pnt = [];
newmrifid.fid.label = {};

D.inv{val}.datareg = struct([]);

% use default selected fiducials
for i = 1:length(meeglbl)
  newmrifid.fid.pnt = [newmrifid.fid.pnt; mrifid.fid.pnt(i, :)];
  newmrifid.fid.label = [newmrifid.fid.label  meeglbl{i}];
end;

M1 = [];
S =[];
S.sourcefid = meegfid;
S.targetfid = newmrifid; 

S.useheadshape = useheadshape; % use headshape points?
ind = 1;

if(do_eeg)
if ~isempty(D.sensors('EEG'))
    if isempty(M1)
        S.template = (D.inv{val}.mesh.template | S.useheadshape);
        M1 = spm_eeg_inv_datareg(S);
    end

    D.inv{val}.datareg(ind).sensors = forwinv_transform_sens(M1, D.sensors('EEG'));
    D.inv{val}.datareg(ind).fid_eeg = forwinv_transform_headshape(M1, S.sourcefid);
    D.inv{val}.datareg(ind).fid_mri = S.targetfid;
    D.inv{val}.datareg(ind).toMNI = D.inv{val}.mesh.Affine;
    D.inv{val}.datareg(ind).fromMNI = inv(D.inv{val}.datareg(ind).toMNI);
    D.inv{val}.datareg(ind).modality = 'EEG';
    
    ind = ind+1;
end
end;

if ~isempty(D.sensors('MEG'))
    if  D.inv{val}.mesh.template
        S.template = 2;
    else
        S.template = 0;
    end

    M1 = spm_eeg_inv_datareg(S);

    D.inv{val}.datareg(ind).sensors = D.sensors('MEG');
    D.inv{val}.datareg(ind).fid_eeg = S.sourcefid;
    D.inv{val}.datareg(ind).fid_mri = ft_transform_headshape(inv(M1), S.targetfid);
    %D.inv{val}.datareg(ind).fid_mri = forwinv_transform_headshape(inv(M1), S.targetfid);
    D.inv{val}.datareg(ind).toMNI = D.inv{val}.mesh.Affine*M1;
    D.inv{val}.datareg(ind).fromMNI = inv(D.inv{val}.datareg(ind).toMNI);
    D.inv{val}.datareg(ind).modality = 'MEG';
end

%spm_eeg_inv_checkdatareg(D);

%D.save;

%%%%%%%%%%%%% forward model

%D = spm_eeg_load(spmfilename);

val=1;
D.inv{val}.forward = struct([]);

ind=1;

if(do_eeg)    
if ~isempty(D.sensors('EEG'))
  D.inv{val}.forward(ind).voltype='EEG BEM';
  ind=ind+1;
end;
end;

if ~isempty(D.sensors('MEG'))
  D.inv{val}.forward(ind).voltype='Single Shell';
  %D.inv{val}.forward(ind).voltype='Single Sphere';

end;

D = spm_eeg_inv_forward(D);

%spm_eeg_inv_checkforward(D, val);

D.save;

if(do_neuromag_grad_baseline_correction)
        D.inv{1}.datareg(1).sensors.tra(1,1:2)
        
        % this bit is needed to correct for error in the fieldtrip code in
        % SPM8
        if(correct_grads)
            correct_factor=1/0.017; % gradiometer coils are 17mm apart
            %correct_factor=1;
            D=correct_planar_grads(spmfilename,['fixed_tra.mat'],correct_factor);
            D.inv{1}.datareg(1).sensors.tra(1,1:2)
        end;
end;
